Research and experimental verification on the mechanisms of cellular senescence in triple-negative breast cancer

Background Triple-negative breast cancer (TNBC) is an aggressive breast cancer subtype with high heterogeneity, poor prognosis, and a low 10-year survival rate of less than 50%. Although cellular senescence displays extensive effects on cancer, the comprehensions of cellular senescence-related characteristics in TNBC patients remains obscure. Method Single-cell RNA sequencing (scRNA-seq) data were analyzed by Seurat package. Scores for cellular senescence-related pathways were computed by single-sample gene set enrichment analysis (ssGSEA). Subsequently, unsupervised consensus clustering was performed for molecular cluster identification. Immune scores of patients in The Cancer Genome Atlas (TCGA) dataset and associated immune cell scores were calculated using Estimation of STromal and Immune cells in MAlignantTumours using Expression data (ESTIMATE) and Microenvironment Cell Populations-counter (MCP-counter), Tumor Immune Estimation Resource (TIMER) and Estimating the Proportion of Immune and Cancer cells (EPIC) methods, respectively. Immunotherapy scores were assessed using TIDE. Furthermore, feature genes were identified by univariate Cox and Least Absolute Shrinkage and Selection Operator (LASSO) regression analyses; these were used to construct a risk model. Additionally, quantitative reverse transcription-polymerase chain reaction (qRT-PCR) and transwell assay were conducted for in vitro validation of hub genes. Result TNBC was classified into three subtypes based on cellular senescence-related pathways as clusters 1, 2, and 3. Specifically, cluster 1 showed the best prognosis, followed by cluster 2 and cluster 3. The levels of gene expression in cluster 2 were the lowest, whereas these were the highest in cluster 3. Moreover, clusters 1 and 3 showed a high degree of immune infiltration. TIDE scores were higher for cluster 3, suggesting that immune escape was more likely in patients with the cluster 3 subtype who were less likely to benefit from immunotherapy. Next, the TNBC risk model was constructed and validated. RT-qPCR revealed that prognostic risk genes (MMP28, ACP5 and KRT6A) were up-regulated while protective genes (CT83) were down-regulated in TNBC cell lines, validating the results of the bioinformatics analysis. Meanwhile, cellular experiments revealed that ACP5 could promote the migration and invasion abilities in two TNBC cell lines. Finally, we evaluated the validity of prognostic models for assessing TME characteristics and TNBC chemotherapy response. Conclusion In conclusion, these findings help to assess the efficacy of targeted therapies in patients with different molecular subtypes, have practical applications for subtype-specific treatment of TNBC patients, and provide information on prognostic factors, as well as guidance for the revelation of the molecular mechanisms by which senescence-associated genes influence TNBC progression.


INTRODUCTION
Triple-negative breast cancer (TNBC), an aggressive breast cancer subtype, lacks the expression of the progesterone receptor, human epidermal growth factor receptor 2 (HER2), and estrogen receptor (Shang & Xu, 2022).As a highly heterogeneous tumor, TNBC shows distinct biological features of aggressiveness, high recurrence rate, and distant metastasis (Gluz et al., 2009;Foulkes, Smith & Reis-Filho, 2010), which can be attributed to germline alterations (BRCA1 mutations), differences in genetic characteristics, epigenetic alterations, DNA repair defects, morphological features, and gene expression profiles (Nik-Zainal et al., 2016;Jiang et al., 2019;Bareche et al., 2018;Pareja & Reis-Filho, 2018).Due to this heterogeneity, large tumors may contain multiple cells with different molecular characteristics and displaying different sensitivity to treatment (Dagogo-Jack & Shaw, 2018), which has been demonstrated to be the main reason for drug resistance or lower efficiency of immunotherapy in breast cancer therapy (Li et al., 2021).Relative to other subtypes, TNBC has a poor prognosis with a low 10-year survival rate of less than 50% (Malorni et al., 2012).Effective therapeutic targets for TNBC are scarce (Huang et al., 2021;Zhu, Yang & Xu, 2022).Therefore, it is particularly crucial to assess the mechanisms underlying the progression and onset of TNBC for immunotherapy and its diagnosis in clinical settings.
The main feature of cellular senescence is a complex secretory phenotype, known as the senescence-associated secretory proteome (SASP), that by facilitate immune-mediated clearance of senescent cells to keep tissue homeostasis and suppresses tumorigenesis, as a powerful tumor-suppressing mechanism (Bartkova et al., 2006;Xue et al., 2007).For example, excessive reactiveoxygenspecies caused by metabolic perturbation could trigger mitogen-activated protein kinase (MAPK)/p38 pathways to facilitate cellular senescence (Borodkina et al., 2016).Researches also uncovered that the MEK/MAPK signaling pathway was involved in inducing senescence in human fibroblasts (Takasugi, Yoshida & Ohtani, 2022).Previous studies suggested that cancer cells have devised various ways to evade this mechanism (Hollstein et al., 1991;Jarrard et al., 1999;Foster et al., 1998).When co-transplanted with completely malignant cells, the growth of senescent cells increases along with the tumor formation rates in the xenografts (Bartholomew, Volonte & Galbiati, 2009;Bhatia et al., 2008;Liu & Hornsby, 2007).The impact of senescent cells on the evolution of precancerous tumor progression has been confirmed (Campisi et al., 2011;Collado, Blasco & Serrano, 2007).Proliferation is induced in several precancerous cells in the presence of senescent fibroblasts, ultimately resulting in tumor formation (Krtolica et al., 2001;Bavik et al., 2006;Castro-Vega et al., 2015).On the one hand, cellular senescence, a fail-safe mechanism, inhibits tumorigenesis in a cell-autonomous environment by preventing harmful cells from further proliferation.On the other hand, senescent cells secrete excess growth factors and cytokines (Coppé et al., 2010), including interleukin 8 (IL-6) and IL-6 (collectively referred to as SASP) (Lasry & Ben-Neriah, 2015), thus establishing a microenvironment of immunosuppression, inflammation, and catabolism to promote chemotherapeutic resistance and tumor growth (Coppé et al., 2008).Intriguingly, within multiple element causing TNBC chemoresistance, therapy-induced senescence is wellaccepted (Chakrabarty et al., 2021).Despite these extensive effects of cellular senescence on cancer, our understanding of cellular senescence-related characteristics of TNBC patients remains rudimentary.In particular, whether these characteristics of cancer patients can be used as biomarkers to guide clinical prognosis and treatment of TNBC remains elusive.
Given that cellular senescence signature is highly correlated with the progression of cancer, an increasing number of studies have investigated it and demonstrated that senescence related signature is significantly associated with the prognosis of tumor patients.For example, in lung adenocarcinoma (LUAD), a comprehensive analysis constructed a senescence-related signature score (SRS), which was elucidated as an independent prognostic factor for LUAD patients and revealed that it was positively associated with cancer-associated fibroblasts, NK cell infiltration, and cytokine release (Lin et al., 2021).Among them, cancer-associated fibroblasts often present a senescence-associated secretory phenotype (SASP) under stress conditions such as inflammatory factors, which promotes tumor formation and chemoresistance through the release of cytokine IL-6, which also reveals the relevant mechanism of cellular senescence-associated cancer progression from the level of regulation of the immune microenvironment (Yasuda et al., 2021;Wu et al., 2017).Another study revealed that SRS in colorectal cancer can be used to assess the genomic mutation status of the tissues, especially the high SRS group showed significant RPN2 mutations (Yue et al., 2021).RPN2 is commonly associated with cell proliferation, migration, and epithelial-mesenchymal transition phenotypes, exacerbates malignant cancer progression through the release of factors such as MMP-9, and inhibits cellular autophagy through STAT3 and NF-κB pathways (Huang et al., 2019;Zhang et al., 2019;Han et al., 2023).Although the clinical value of cellular senescence-related features has not been explored in TNBC, the above reports still inform our study, and we hypothesize that markers associated with senescence features and prognosis in TNBC could influence cancer progression and prognosis in terms of tumor infiltration microenvironment homeostasis, cytokine release, and genomic mutations.
Based on scRNA-seq and bulk RNA-seq data, we constructed a senescence classifier for TNBC based on cellular senescence-related pathways.Significant differences in prognosis, tumor microenvironment (TME), genome mutations and enriched metabolic pathways were detected from the three subtypes.Next, the heterogeneity of malignant cells among senescent and normal cells was further analyzed in TNBC.A cellular senescence-related risk model for TNBC was constructed.In vitro experiments were also conducted to confirm the expression levels of hub genes selected for the risk model and the functions of certain gene on TNBC progression.Additionally, we also verified the reliability of the signature in assessing the regulatory role of the TNBC immune microenvironment and the patient's response to immunotherapy.Finally, the possibility of cellular senescence as a clinicopathological feature to improve the prognosis and survival prediction of patients with TNBC was explored.Collectively, our study provides a new direction to the diagnosis and treatment of patients with TNBC.

Data downloading and pre-processing
Single-cell sequencing dataset, GSE176078, comprising nine samples was extracted from the NCBI data Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/).The statistical power (Therneau, Hart & Kocher, 2023) of this experimental design, calculated in RNASeqPower is 1.0.Single-cell data were subjected to filtering, quality control, clustering, and dimensionality reduction using the ''Seurat'' (Gribov et al., 2010) package.Microarray data of GSE58812 was extracted from GEO, and probes were transformed to ''symbols'' using the annotation file.In total, 107 qualified tumor samples and 16,416 genes were obtained after excluding tumor samples and normal tissues without information on clinical follow-up and overall survival (OS).
Clinical phenotypic data of TNBC were extracted from The Cancer Genome Atlas (TCGA) database.Samples without information on survival time and statuses were removed, and all those with patient survival time >0 days were retained.Subsequently, expression profile data were downloaded from TCGA.Finally, 113 para-cancerous and 113 tumor samples were obtained.The statistical power of this experimental design, calculated in RNASeqPower is 1.0.Furthermore, copy number variations (CNVs) on the Masked Copy Number Segment type of TNBC were obtained from TCGA database and integrated using ''gistic2'' (Mermel et al., 2011) software.Finally, the single nucleotide variants (SNVs) in TCGA-TNBC dataset processed by ''mutect2'' software, were obtained from TCGA database.Cellular senescence-related pathways were extracted from the Molecular Signature Database (MSigDB: https://www.gsea-msigdb.org/gsea/index.jsp).Sangerbox platform (http://vip.sangerbox.com/) was introduced to assist bioinformatics analysis in this study (Shen et al., 2022).

Single-cell clustering and dimensionality reduction
First, 38,582 cells were obtained after filtering the single-cell RNA-seq data matrix according to the cell criterion (<250 transcripts/cell) and gene criterion (<3 cells/gene).The ''PercentageFeatureSet'' function was then used to determine the percentages of mitochondria and rRNAs.Ultimately, 38,007 cells that met the following inclusion criteria were obtained: One cell could express between 100 and 6,000 genes, the mitochondrial content was less than 25%, and each cell had a unique molecular identifier (Yasuda et al., 2021) of at least 100.
Variance stabilization transformation (VST) was used to identify variable characteristics after the aforementioned data were log-normalized.High-variance genes were then discovered by employing the ''FindVariableFeatures'' tool.Further, the ''FindIntegrationAnchors'' function was utilized to remove batch effects using the canonical correspondence analysis (CCA) method for 10 samples.The data were integrated using the ''IntegrateData'' function.The anchor (dim = 30) was determined by principal component analysis (PCA) and dimensionality reduction after normalization using the ''ScaleData'' function for all genes.The cells were clustered using ''FindNeighbors'' and ''FindClusters'' functions with resolution = 0.1; this yielded 11 subpopulations.Finally, the TSNE dimensionality reduction analysis was performed on 38,582 cells using the ''RunTSNE'' function.
With threshold parameters of logfc = 0.5 (fold-change), Minpct = 0.35 (minimum percentage of differential gene expression), and adjusted P0.05, marker genes of subpopulations were screened using the ''FindAllMarkers'' tool.Subsequently, changes in cellular CNVs in single-cell data were predicted using the ''copycat'' (Gao et al., 2021) package to distinguish tumor cells.Although cancerous and normal tissues were differentiated during sampling, cancerous tissues might contain normal cells.Therefore, a distinction between the two was necessary.

Gene set enrichment analysis (GSEA) and annotation
Cellular senescence-related pathways were obtained from the GSEA website (https: //www.gsea-msigdb.org/gsea/index.jsp).We separately calculated the corresponding scores of non-malignant and malignant cells in cellular senescence-related pathways by ssGSEA using the ''GSVA'' package and normalized the enrichment scores for each pathway by z-score.The scores of each sample of normal and tumor tissues in the bulk RNA-seq dataset were analyzed by ssGSEA for these senescence-related pathways (Subramanian et al., 2005).Finally, the significance of each senescence-related pathway in cancer and para-cancerous tissues was assessed using ''wilcox.test''.
The genes related to the G1/S phase (p15-Cell cycle G1/S, MDM2-p21-Cell cycle G1/S, and p27-Cell cycle G1/S; 27 genes in total) were extracted from the Kyoto Encyclopedia of Genes and Genomes (KEGG) official website.Subsequently, the scores for G1/S genes were calculated for each sample in TCGA dataset by ''ssGSEA''.Genes in the HALLMARK_G2M_CHECKPOINT pathway were downloaded from MSigDB via GSEA and the scores of G2 checkpoints were calculated for each sample in TCGA dataset by ''ss-GSEA''.We downloaded REACTOME_TELOMERE_EXTENSION_BY_TELOMERASE using GSEA and the primary role of this pathway was ''Telomere Extension By Telomerase''.Finally, GSEA-derived enrichment scores in TCGA dataset for epithelialmesenchymal transition (EMT) were computed by ssGSEA for 200 genes in the HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION pathway for each sample (Liberzon et al., 2015).

Univariate Cox and LASSO regression analyses
Univariate Cox analysis was conducted using the ''coxph'' function in the survival package of R to screen the significant (p < 0.05) prognosis-related genes.LASSO-Cox regression analysis was conducted for prognosis-related genes using the ''glmnet'' package in R (Friedman, Hastie & Tibshirani, 2010).The changing trajectory of each independent variable was analyzed.Increased lambda was positively related to increased number of independent variable coefficients tending to zero.Finally, by 10-fold cross-validation, a model was constructed; for each lambda, the confidence intervals were calculated.
The risk-related prognostic score (Riskscore) was calculated for each sample according to the equation for the sample risk score: Riskscore = βi×Expi, where Expi denotes the level of gene expression of the corresponding gene signature and β is the corresponding Cox regression coefficient.ROC analysis and z-score for RiskScore were computed using the ''timeROC'' package in R. Samples with RiskScore lower than zero comprised the low-risk group, while those with RiskScore greater than zero after z-score normalization were classified into the high-risk group.Finally, the Kaplan-Meier curves were plotted.

Consensus clustering
Consensus clustering was conducted using the ''ConsensusClusterPlus'' package (Wilkerson & Hayes, 2010).In total, 500 bootstraps were conducted using the ''hc'' algorithm with ''canberra'' as the metric distance; in the training set, 80% of the patients were in each bootstrap process.After the number of clusters was selected to between 2 and 10, ''ConsensusClusterPlus'' in the TCGA clustered the 113 TNBC samples dataset.The optimal classification was confirmed based on the consensus matrix and cumulative distribution function (CDF).

Mutation analysis
Using hg38 as the reference genome, the CNV-related TCGA-TNBC results were integrated by ''gistic2 software at a confidence level of 0.9.The downloaded SNV data from TCGA were analyzed using the ''maftools'' (Mayakonda et al., 2018) package.

Immune cell scoring and immunotherapy
The scores of the relevant immune cells were calculated using the ''MCPcounter.estimate''function of the ''MCPcounter'' package (Becht et al., 2016), Tumor Immune Estimation Resource (TIMER) (Li et al., 2020) and Estimating the Proportion of Immune and Cancer cells (EPIC) (Zhou et al., 2020) methods.Significant differences were determined using ''kruskal.test''.Immune cell infiltration was assessed using the ''Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression'' (ESTIMATE) (Yoshihara et al., 2013) which provides information on tumor purity, scores of immune cell infiltration, and stromal cell levels in tumor tissues.Subsequently, the TIDE (Jiang et al., 2018) software was used to assess the putative immunotherapeutic efficacy in the defined molecular subtypes (http://tide.dfci.harvard.edu/).Higher TIDE score indicates less immunotherapy benefit as it suggests a greater possibility of immune escape.The correlation and significance of the risk score with the immune cell score were separately calculated using the ''Hmisc'' package's ''rcorr'' function based on Pearson's method.

Quantitative reverse transcription-polymerase chain reaction
Using TRIzol reagent (Thermo Fisher, Waltham, MA, USA), total RNA was collected from MDA-MB-468, MDA-MB-231 and MCF10A cell lines with 260/280 the value was between 1.8-2.0.High-Capacity cDNA Reverse Transcription Kit (4368814, ThermoFisher, Waltham, MA, USA) was conducted for reverse transcription.Using a LightCycler 480 PCR System (Roche, Indianapolis, IN, USA), the RNA from each sample (2 µg) was subjected to qRT-PCR with FastStart Universal SYBR ® Green Master (Roche, Indianapolis, IN, USA).The cDNA was a template with a reaction volume of 20 µl (0.5 µl of forward and reverse primers, 2 µl of cDNA template, 10 µl of PCR mixture, and appropriate volume of water).For the PCR reactions, cycling began with an initial DNA denaturation phase for 30 s (s) at 95 • C, then a total of 45 cycles were run for 15 s at 94 • C, for 30 s at 56 • C, and for 20 s at 72 • C. Each sample was performed for three separate analyses.Data from the threshold cycle (CT) were normalized to the level of GAPDH by 2 − CT .The mRNA expression was compared in the normal and control tissues.Sequences of primer pairs targeting (Shanghai Gemma Gene Co., LTD, Shanghai, China) the genes are listed as follows:

Transwell assay
Migration and invasion of MDA-MB-468 and MDA-MB-231cell lines were detected.Cell (5 × 10 4 ) inoculation onto chambers coated (for invasion) or uncoated with Matrigel (BD Biosciences, USA) (for migration) was performed.The upper layer was added with serum-free medium, while the lower layer was added with complete DMEM medium.After incubation for 24 h, migrated or invaded cells were fixed using 4% paraformaldehyde and dyed by 0.1% crystalline violet.Cell number was counted using a light microscope.The experiment was conducted in our own lab.

Statistical analysis
All statistical analysis was done with R program (v4.2.0).The Wilcoxon test evaluated twogroup difference.Differences among three groups were analyzed by the Kruskal-Wallis test.
Cox analysis and survival analysis both used the log-rank test.For Wet experiments, the number of each group was 3 independent experiments.Assay carried out by investigator's lab.Consumables in experiments we use are RNAase and DNAase free.

Nine annotated cell types displayed different senescence characteristics
A cellular senescence-related classifier was constructed for TNBC.The prognoses of patients with TNBC were verified using bulk RNA-seq and scRNA-seq data.The flow chart of the study design is presented in Fig. S1.First, 38,582 cells were obtained by screening the data from nine samples (Table S1, Fig. 1A); these were clustered and 11 subpopulations were obtained (Fig. 1B).Subsequently, these cells were subjected to TSNE dimensionality reduction analysis, and the 11 subpopulations were annotated based on some classical markers of immune cells (Fig. 1C and Table S2).Among them, subpopulation C1 comprised CD8 T cells expressing CD8A, CD3D, GZMA, and CD8B; C6 comprised CD4 T cells expressing CD4 and CD3D; C2 and C9 comprised macrophages expressing CD163 and CD68; C3 and C10 comprised monocytes expressing S100A8; C5 comprised B cells expressing CD19, CD79A, and MS4A1; C4 comprised plasma cells expressing CD79A and JSRP1; C0 comprised epithelial cells expressing EPCAM; C3 and C8 contained fibroblasts expressing ACTA2, FAP, PDGFRB, and NOTCH3, and C7 contained endothelial cells expressing PECAM1.We also screened marker genes among these subpopulations and analyzed their corresponding expression (Fig. 1D).For instance, CD24 and KRT19 were highly expressed in epithelial cells in C0 and were the marker genes of C0.
The proportion of these nine subpopulations in each sample was further analyzed (Fig. 1E).Subsequently, changes in CNVs in single cells were predicted to distinguish normal from tumor cells in each sample (Fig. 1F).In total, 7,709 cancer cells and 30,298 normal cells were identified.Subsequently, the proportion of malignant and non-malignant cells in each sample was calculated (Fig. 1G).The cells in samples CID3946, CID44041, and CID4465 were all non-malignant, while the other samples comprised both malignant and non-malignant cells.
The senescence characteristics of single cells were further assessed.Cellular senescencerelated pathway scores were higher in malignant relative to the non-malignant cells (Fig. 2).

Validation of abnormal cellular senescence based on bulk RNA-seq data
Cellular senescence-related pathways were highly expressed in malignant cells at the single cell level relative to the non-malignant cells, with aberrant expression.Therefore, their expression in tumor and normal tissue samples was further analyzed using a bulk RNA-seq dataset.GSEA results suggested that REACTOME_CELLULAR_SENESCENCE, and P53_SIGNALING_PATHWAY were significantly enriched in tumor tissues in TCGA dataset (Fig. 3A).ssGSEA of the bulk RNA-seq dataset (Fig. 3B) revealed that the enrichment scores of REACTOME_DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE, P53_SIGNALING_PATHWAY, REACTOME_CELLULAR_SENESCENCE, REAC-TOME_DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE were higher in some cancer tissues relative to the para-cancerous tissues.Therefore, these three cellular senescence-related pathways may exhibit regulatory complexity in cancerous tissues.

Construction of cellular senescence-related subtypes
Three cellular senescence-related pathways, including REACTOME_CELLULAR_ SENESCENCE, and P53_SIGNALING_PATHWAY, REACTOME_DNA_DAMAGE_ TELOMERE_STRESS_INDUCED_SENESCENCE, were significantly enriched in TNBC tissues.Therefore, we further analyzed the genes enriched in these three pathways.The above pathways comprised 253 genes, of which 186 were present in TCGA dataset.Subsequently, a univariate analysis of these 186 genes yielded seven prognosis-related genes, including ETS2, SERPINE1, FOS, SHISA5, IL1A, TP53AIP1, and IGFBP7 (Fig. 4A).Mutations in these seven genes in TNBC were further examined.Among SNVs, ETS2 showed the highest mutation frequency (primarily mutation type: missense) (Fig. 4B).Among CNVs, IGFBP7 had the highest frequency of acquired CNVs (CNV_gain), whereas SHISA5 had the highest frequency of deletion mutations (CNV_loss) (Fig. 4C).The levels of these seven genes in normal versus tumor tissues were analyzed by the ''wilcox.test''(Fig. 4D).The expressions of SHISA5, SERPINE1, and IL1A were markedly high in tumor tissues, whereas those of TP53AIP1, ETS2, and FOS were high in the normal tissues.
Consensus clustering of 113 TNBC samples in TCGA dataset was performed based on these seven key prognosis-related genes.Stable clustering results were achieved at k = 3 (Figs.4E-4F).Therefore, three subtypes (clusters) were obtained (Fig. 4G).Significant differences in the prognostic characteristics among the subtypes were found (Fig. 4H).Overall, cluster 1 showed the best prognosis, followed by clusters 2 and 3.These seven genes showed the highest expression in cluster 3 (Fig. 4I).

Differential analysis of cellular senescence-related subtypes
The expression of genes in these three senescence-related pathways among the three clusters was compared by ''kruskal.test''(Fig. S2).Among the three pathways, genes in cluster 2 showed the lowest expression, while those in cluster 3 showed the highest expression.
The distribution of clinical characteristics across molecular subtypes was compared by chi-square test to assess the differences among the three subtypes in TCGA.The survival statuses of these patients differed significantly among the three subtypes.In particular, the proportion of dead patients was high in cluster 3. Differences in mutations among cellular senescence-related subtypes were analyzed.TCGA data suggested some differences in CNVs among the three subtypes (Figs.5A and  5C).For instance, cluster 1 had the lowest degree of gene gain and lose.While cluster 3 exhibited the highest degree of CNV gain.Additionally, clusters 2 and 3 had similar CNV loss on chromosomes 2, 3, 4 and 5. Furthermore, TP53, TTN, MUC16, SYNE1, and FAT3 had the highest SNV frequency (Fig. 5B).

Biological characteristics of cellular senescence-related subtypes
Cancer cells can induce cellular senescence through the inhibition of the cell cycle, and cell cycle protein-dependent kinases (CDK inhibitory proteins), such as InK4a and p21, are upregulated in senescent cells, which causes cell cycle arrest (Wang, Lankhorst & Bernards, 2022).A previous study Cuzick et al. (2011)  progression (CCP).We computed the CCP scores of each sample in TCGA dataset by ssGSEA.clusters 1 and 2 had slightly higher CCP scores relative to cluster 3 (wilcox.test)(Fig. 6A).No significant differences were observed among the three subtypes in the G1/S phase (wilcox.test)(Fig. 6B), whereas the scores for G2 checkpoints were slightly lower in cluster 3 relative to clusters 1 and 2 (wilcox.test)(Fig. 6C).Prognostic analysis suggested that patients in cluster 3 showed a worse prognosis.Therefore, the cell cycle may be one of the factors influencing cellular senescence.Other mechanisms may also regulate cellular senescence along with the cell cycle.
Similarly, inhibition of telomerase induces cellular senescence (Wang, Lankhorst & Bernards, 2022).Cancer cells usually circumvent telomere attrition by activating telomerase activity.Telomere elongation scores for telomerase in cluster 3 were lower than those in Full-size DOI: 10.7717/peerj.16935/fig-5 clusters 1 and 2, with the latter showing the highest score.Therefore, other mechanisms may co-regulate cellular senescence with telomerase activity (wilcox.test)(Fig. 6D).However, in addition to conveying the message ''please kill me'', the factors secreted by senescent cells also affect surrounding cells (Yu et al., 2021).This effect can promote tumor migration and metastasis by promoting EMT.Senescent tumor cells can promote the production of blood and lymphatic vessels by recruiting specialized macrophages that provide oxygen and nutrients necessary for growth to other tumor cells, thus promoting tumor growth and metastasis.We calculated the EMT scores.Cluster 2 had a lower EMT score than clusters 1 and 3, indicating that patients in cluster 3 were more likely to experience metastasis (wilcox.test)(Fig. 6E).
Based on genes in the HALLMARK_HYPOXIA pathway, the hypoxia scores of the samples were computed by ssGSEA.Angiogenesis scores of samples were also analyzed based on 24 genes from the literature (Masiero et al., 2013) (wilcox.test)(Figs.6F-6G).Enrichment scores for 10 pathways associated with tumors in each sample in TCGA dataset were computed by the ssGSEA method described previously (Sanchez-Vega et al., 2018).Only two of the ten pathways (NOTCH and TGF-beta signaling) exhibited significant differences (wilcox.test).
Finally, genes in the HALLMARK_INFLAMMATORY_RESPONSE pathway from GSEA-based MSigDB were analyzed for inflammation-related scores by ssGSEA.Cluster 2 had a significantly lower inflammation-related score relative to clusters 1 and 3 (wilcox.test)(Fig. 6H).

Immune characteristics of cellular senescence-related subtypes
The relationship between cellular senescence-related subtypes and immunity was analyzed to investigate the immune characteristics of these subtypes.The immune scores of patients in TCGA dataset were calculated by ''ESTIMATE''.Significant differences were observed among the three subtypes, with a high degree of immune infiltration in clusters 1 and 3 (kruskal.test)(Fig. 7A).Subsequently, the scores for the relevant cells were calculated.Significant differences were noted for all cells except CD8 T cells and cytotoxic lymphocytes (kruskal.test)(Fig. 7B).Similarly, T cells CD4, also displayed distinct difference among three groups in TIMER and EPIC analysis (Figs. 7C-7D).Although the results were obtained from different algorithms, these results all together suggested higher immune infiltration status in clusters 3. The expression of immune checkpoint genes, including NRP1, CD200, BTLA, LAIR1, and TNFRSF14 differed significantly among the three subtypes (kruskal.test)(Fig. 7E).Finally, the putative clinical efficacy of immunotherapy was assessed among these molecular subtypes.High TIDE scores in cluster 3 in TCGA cohort suggested that patients of this subtype were less likely to benefit from immunotherapy owing to a higher likelihood of immune escape (wilcox.test)(Fig. 7F).

Identification of prognosis-related genes and construction of the risk model
Differential analysis was separately performed between cluster 1 and non-cluster 1; cluster 2 and non-cluster 2, and cluster 3 and non-cluster 3. Finally, 2 down-regulated and 34 up-regulated DEGs were screened between cluster 1 and non-cluster 1; 293 down-regulated and 9 up-regulated DEGs between cluster 2 and non-cluster 2, and 26 down-regulated and 218 up-regulated DEGs between cluster 3 and non-cluster 3. Thus, a total of 391 DEGs were screened.
The four genes' expression levels were utilized to determine the risk score for each sample using the TCGA dataset as the training set.Then, using ROC analysis, the RiskScore prognostic classification was carried out (Fig. 8A).The prognostic prediction classification efficiency was analyzed separately for 1-5 years with obtained AUC values 0.91, 0.93, 0.74, 0.83 and 0.84, respectively , and the AUC values reached 0.7, indicating good predictive performance.Samples with RiskScore greater than zero after z-score normalization were classified into the high-risk group, whereas those with Riskscore lesser than zero comprised the low-risk group, with a highly significant (p < 0.05) difference between them (Fig. 8B).
To validate the model's robustness, the GSE58812 dataset was used for verification.Specifically, the risk model was constructed using these four genes and ROC analysis for prognostic classification was performed based on the RiskScore (Fig. 8C).The prognostic predictive classification efficiency was analyzed for 1-5 years with obtained AUC values 0.98, 0.71, 0.62, 0.66 and 0.65, respectively.All AUC values were over 0.6, displaying good predictive ability.Furthermore, a significant difference in prognosis was found between the risk groups (p < 0.05) (Fig. 8D).

Validation of association between prognostic model signature and TNBC malignant phenotype
In this study, molecular assays revealed that the expression of MMP28, ACP5 and KRT6A was increased in MDA-MB-468 and MDA-MB-231 cell lines.In contrast, the expression of CT83 was decreased in MDA-MB-468 and MDA-MB-231 cell lines (Figs. 9A-9D).This result is consistent with previous analyses that high expression of prognostic risk factors and low expression of protective factors associated with higher RiskScore.Given that ACP5 is a key gene influencing high RiskScore, we inhibited the expression of ACP5 in MDA-MB-468 and MDA-MB-231 cell lines, and by transwell assay, we could observe the ability of ACP5 to promote the proliferation of triple-negative breast cancer cell lines.After we inhibited the expression of ACP5 in MDA-MB-468 and MDA-MB-231 cell lines, the invasion and migration ability of the cell lines were decreased (Figs.9E-9G).This study not only confirms the existence of ACP5 as a prognostic risk factor, but also demonstrates that this gene regulates the TNBC malignant phenotype at the cellular level, thereby influencing TNBC malignant progression.

Immunological characteristics of RiskScore subgroups and sensitivity analysis based on drug IC50
To elucidate differences in the immune microenvironment of patients between the risk groups, differences in the relative cell abundances predicted by ''MCPcounter'' were compared (wilcox.test)(Fig. 10A).T cells, cytotoxic lymphocytes, and monocytic lineage differed significantly between the risk groups.Immune cell infiltration was also assessed (wilcox.test)(Fig. 10B).''ImmuneScore'' was higher in the high-risk group relative to the low-risk group, suggesting that patients in the former had high immune cell infiltration.Additionally, TIMER and EPIC analysis also revealed higher immune infiltration status in high risk group (Figs.10C-10D).Subsequently, differences in response to immunotherapy between the risk groups were analyzed in TCGA cohort.First, differences in immune checkpoint expression were compared.Most immune checkpoint genes, including TNFRSF14, LAIR1, and CD244 showed differential expression between the risk groups, with high expression in the high-risk group (wilcox.test)(Fig. 10E).
A significant positive correlation between the Riskscore and T cells, cytotoxic lymphocytes, and monocytic lineage (Fig. 10F) was observed.Risk scores also correlated positively with ''ImmuneScore'', ''StromalScore'', and ''ESTIMATEScore'' (Fig. 10G).Drug sensitivity prediction revealed that low-risk group showed higher sensitivity to six traditional drugs (MG-132, Sorafenib, A77004, Bortezomib, Shikonin and AZ628) by wilcox.test(p < 0.05, Fig. 10H).Univariate regression analysis based on RiskScore and other clinicopathological characteristics suggested that the latter was the most significant prognostic factor (pvalue = 0.007, Fig. 11A) and remained so (p-value = 0.003, Fig. 11B).A nomogram was constructed to quantify the survival probability of patients and risk assessment by combining RiskScore and other clinicopathological characteristics (Fig. 11C).The results confirmed the greatest impact of RiskScore on survival prediction.Subsequently, the model's prediction accuracy was assessed using a calibration curve (Fig. 11D).The findings indicated that the calibration curves for 1-, 3-, and 5-year predictions nearly overlapped with the standard curves, pointing to the nomogram's strong predictive ability.Furthermore, decision curve analysis (DCA) was used to evaluate the model's dependability.DCA results highlighted that both RiskScore and nomogram had significantly more benefits than the extreme curves with strong survival predictive power relative to other clinicopathological features, especially for the 5-year survival (Figs.11E-11F).

DISCUSSION
TNBC is a highly heterogeneous type of cancer (Bianchini et al., 2016).Recent studies have classified TNBC based on different molecular or clinical subtypes and identified possible therapeutic targets, primarily based on transcriptomic subtypes (Metzger-Filho et al., 2012;Garrido-Castro, Lin & Polyak, 2019;Hu et al., 2022).Cellular senescence promotes tumor regression through cell-non-autonomous and -autonomous mechanisms.Drugs inducing cancer cell senescence and modulating SASPs provide new targets for tumor therapy (Calcinotto et al., 2019).However, no cellular senescence-based typing of TNBC has been established to date.In this study, we proposed cellular senescence molecular typing for the first time, and identified seven genes related to TNBC prognosis through analyses and established three molecular subtypes.Subsequently, DEGs among three molecular subtypes was adopted for the construction of risk model.Finally, a nonogram for clinical decision-making was designed.
Seven cellular senescence related genes exhibited different biological functions in TNBC or other cancer development.ETS2 is a key ETS transcriptional family member implicated in invasion and metastasis in TNBC (Kollareddy & Martinez, 2021).A screening experiment conducted against CDK10 confirmed the ETS2 transcription factor as an interacting protein (Guen et al., 2017).This finding indicated that ETS2 may also participate in cell cycle process.Serpin family E member 1 (SERPINE1) is responsible for encoding plasminogen activator inhibitor 1, which is a primary inhibitor of tissue plasminogen activator (Li et al., 2018).SERPINE1 has been detected in various cancer and involved in cancer invasion, migration, and angiogenesis (Seker et al., 2019;Yang, Ma & Zhu, 2019).Also, SERPINE1, a paclitaxel-resistant oncogene, is a putative target for TNBC treatment (Zhang, Lei & Jing, 2020b).FOS, a key regulator of TNBC cell proliferation and viability, is related to a poor patient prognosis (Zhang, Lei & Jing, 2020a).IL1A, a pro-inflammatory cytokine, enhances cell growth and invasiveness of TNBC cells.Abnormal IL1A induction correlates with a poor patient prognosis in TNBC (You et al., 2021).IGFBP7 was reported to suppress breast tumor growth through the induction of apoptosis and senescence pathways (Benatar et al., 2012).Overall, these findings disclosed the potential functions or the prognostic ability of selected hub genes and supported the reliability of our results.
TNBC is the most aggressive and heterogeneous subtype of breast cancer, and consists of several different microenvironmental phenotypes (Xiao et al., 2019;Gruosso et al., 2019;Bareche et al., 2020).Previous studies have classified several molecular subtypes for TNBC.For instance, Lehmann and co-workers carried out clustering analysis and categorized TNBC patients into six subtypes, called basal-like 1, basal-like 2, immunomodulatory, luminal androgen receptor, mesenchymal, and mesenchymal stem-like subclusters (Lehmann et al., 2011).Next, a transcriptome-wide analysis of Chinese TNBC samples clustered TNBC patients into four subtypes -BLIS, IM, LAR, and MES (Tong et al., 2023).In this research, we classified TNBC patients into only three clusters based on senescence related characteristics, and our immunological analysis highlighted significant differences in TME among the three subtypes, which could support TNBC cellular senescence-based therapy and support clinical decision-making.Changes in the vital immune cells that make up the stroma around the tumor are correlated with tumor growth.TME metabolism changes are the result of the aberrant metabolic status of tumor cells (Hinshaw & Shevde, 2019).The immune score produced by the ESTIMATE algorithm rates the immunological make-up and speed of tumor sample responses.Strong correlations exist between prognosis and tumor purity, or the proportion of malignant cells in the tumor tissue (Aran, Sirota & Butte, 2015).Here, we used ''ESTIMATE'' to identify three subtypes in TME.Significant differences were found among the three subtypes, with patients in clusters 1 and 3 showing a high degree of immune infiltration.Actually, senescent cells secrete various factors which make uo senescence-relevant secretory phenotype, including pro-inflammatory factors (cytokines, chemokines, micro-RNAs) to induce inflammatory state (Frasca et al., 2021).Therefore, in this work, TME-related cellular senescence characteristics are of interest in determining treatment strategies for TNBC and could be potential targets for individualized therapy.
We built a risk model based on genes in cellular senescence-related pathways and calculated risk scores based on four significant genes, including MMP28, CT83, ACP5, and KRT6A, and validated the relative expression levels of the four genes in TNBC cells.A previous study defined a subset of TNBC with poor prognosis based on the basal marker, KRT16, the stem cell marker, WNT11, and the EMT marker, MMP28 (Yu et al., 2013), corroborating the rationality and validity of our gene selection approach.We constructed a risk model incorporating multiple genes, wherein MMP28 was a risk factor.EMT scores were calculated.Patients in cluster 2 had lower EMT scores than those in clusters 1 and 3, suggesting that the former were more likely to experience metastasis.CT83 is a gene specific to TNBC and its hypermethylation is oncogenic in breast cancer (Chen et al., 2021).Overexpression of four genes, including ART3, FABP7, TTYH1, and CT83 correlated positively with patients' life expectancy rates in TNBC (p < 0.05) (Zhong et al., 2020).This is supported by our study revealing the down-regulated expression of this gene in the high-risk group.A model comprising CT83 and FABP7 constructed using Cox regression, Kaplan-Meier, and ROC analyses correlated significantly with OS in TNBC.ACP5 has been shown to correlate with poor patient prognosis in cancers such as gastric and rectal cancers, and to correlate with malignant tumor progression by affecting malignant phenotypes such as cell proliferation and invasion (Bian et al., 2019;An et al., 2021).In breast cancer-related studies, the TRAP/ACP5/uteroferrin/purple acid phosphatase/PP5 signaling axis can act as a driver mediating breast cancer invasion and mediate cancer malignant progression at the cellular level (Krumpel et al., 2015).This informs the present study, in which we reveal that this gene is up-regulated for expression in TNBC cells and promotes cell migration and invasive ability, providing a reference for subsequent work to uncover the mechanisms by which this gene affects TNBC progression at the cellular level.Ultimately, we identified potential biomarkers implicated in the onset and progression of TNBC.These genes may provide novel insights into the tumorigenesis of TNBC and serve as independent prognostic factors for TNBC.Taken together, these findings demonstrate the validity of the risk model developed in this study.
Some limitations exist in this study that should be stated prior to outlining the conclusions.Firstly, we could not demonstrate the difference and role of these subtypes in TNBC progression because of a lack of progression data such as tumor stage of the patients.Second, we only verified the role of the key signature gene associated with RiskScore in the regulation of the malignant phenotype of TNBC cells at the cellular level, but we did not verify the molecular mechanisms of the signature gene, such as influencing cytokine secretion to regulate cellular senescence, to corroborate the results of the bioinformatic analyses.Finally, rather than evaluating our cohort, the data were taken from available databases.As a result, there was still little evidentiary impact.Prospective studies are thus needed to validate the therapeutic and prognostic value of the established subtypes in TNBC.

CONCLUSION
In summary, we elucidated the abnormalities among cellular senescence-related pathways in TME of malignant cells of TNBC.Senescence-related subtypes and a risk model were constructed based on the genes in these pathways.Also, we validated the oncogenic gene ACP5's ability to migrate and invade in TNBC cells.Additionally, a nomogrom was designed for future clinical decision-making.All these findings provide a basis for the revelation of the molecular mechanisms by which senescence-associated genes influence TNBC progression and directions for TNBC therapy and prognosis.

Figure 1
Figure 1 TSNE dimensionality-reduction results.(A) TNSE plots for the distribution of nine samples; (B) TSNE plot for 11 immune cell subpopulations; (C) TSNE plot for cell distribution after annotation; (D) dot plots of the top five marker genes' expression in the annotated subpopulations; (E) proportions and cell count of the annotated subpopulations in each sample; (F) distribution of malignant and nonmalignant cells predicted by ''copykat''; (G) proportion of malignant and non-malignant cells in each sample.Full-size DOI: 10.7717/peerj.16935/fig-1

Figure 3
Figure 3 Enrichment analysis for bulk RNA-seq data.(A) GSEA results in TCGA dataset; (B) Heat map of the distribution of cellular senescence-related pathway ssGSEA scores in cancerous and para-cancerous tissues in TCGA dataset.Full-size DOI: 10.7717/peerj.16935/fig-3

32 Figure 4
Figure 4 Analysis of cellular senescence-related prognostic genes.(A) Forest plot of prognosis-related genes; (B) waterfall plot of SNVs in the seven prognosis-related genes; (C) Percentage plot of CNVs in the seven prognosis-related genes; (D) The expression of the seven prognosis-related genes expressed in tumor and normal tissues; (E) CDF curve of TCGA cohort; (F) CDF-delta area curves for TCGA cohort.The Delta area curve for consensus clustering indicates the relative change in the area under the CDF curve for each category number, k, relative to k-1.The horizontal axis represents the category number, k, and the vertical axis represents the relative change in area under the CDF curve; (G) Heat map of sample clustering at consensus k = 3; (H) KM curve of the relationship between the prognoses of the three subtypes in TCGA; (I) Heat map of the expression of seven prognosis-related genes among three subtypes in TCGA dataset.Full-size DOI: 10.7717/peerj.16935/fig-4

Figure 5
Figure 5 Differential analysis of cellular senescence-related subtype mutations.(A) Comparison of CNVs among subtypes; (B) Waterfall plots of the top 10 genes with the highest SNVs among subtypes; (C) Peak plots of genes amplified (red) and missing (blue) among the three subtypes (G scores (top) and qvalues (bottom) relative to the entire graph for amplification of all markers on the analyzed region).Full-size DOI: 10.7717/peerj.16935/fig-5

Figure 8
Figure 8 Construction and validation of the risk model.(A) ROC curves for the risk model constructed using four genes from TCGA dataset analysis; (B) KM curves for the risk model based on TCGA dataset; (C) ROC curves for the risk model constructed using four genes for the GSE58812 dataset; (D) KM curves for the risk model based on the GSE58812 dataset.Full-size DOI: 10.7717/peerj.16935/fig-8